#include <emmintrin.h>
#include <omp.h>

int min(int a, int b);
void single_edgecase(int m, int jboun, int kboun, int iboun, float *A, float *A2, float *C);

int min(int a, int b)
{
  if(a<b)
     return a;
  return b;
}

void sgemm( int m, int n, float *A, float *C )
{

  __m128 j1,j2,j3,j4,c1,i_col1, i_col2, i_col3, i_col4, i_col5;

  int i4, i8, i12;

  float *Ctemp;  
  float *Ctemp2;
  float *Ctemp3;
  float *Ctemp4;
  #pragma omp parallel for private(j1,j2,j3,j4,c1,i_col1, i_col2, i_col3, i_col4, i_col5, i4, i8, i12, Ctemp, Ctemp2, Ctemp3, Ctemp4  )
  
 for(int jj=0;jj<m;jj+=64) {
  for(int kk=0; kk<n;kk+=64) {
   for(int ii=0;ii<m;ii+=64) {
  int iboun = min(64, (m-ii));
  int jboun = min(64, (m-jj));
  int kboun = min(64, (n-kk));
  if(iboun==64&&jboun==64&&kboun==64) {
  for( int k = kk; k < kk+64; k++ )
  {
    
    for( int j = jj; j < jj+64; j+=4 )
    {
      int stored = j+1, stored2=j+2, stored3=j+3;
      j1 = _mm_load1_ps(A+(j+k*m));
      j2 = _mm_load1_ps(A+((stored)+k*m));
      j3 = _mm_load1_ps(A+((stored2)+k*m));
      j4 = _mm_load1_ps(A+((stored3)+k*m));

      Ctemp = C+(j*m);
      Ctemp2 = C+(stored*m);
      Ctemp3 = C+(stored2*m);
      Ctemp4 = C+(stored3*m);
      for( int i = ii; i < ii+64; i+=16 )
      {
        //C[i+j*m] += A[i+k*m] * A[j+k*m];
	i4=i+4;
	i8=i+8;
	i12=i+12;
        //C[i+j*m] += A[i+k*m] * A[j+k*m];
        c1 = _mm_loadu_ps(Ctemp+i);
        i_col1 = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j1));
        _mm_storeu_ps(Ctemp+i,c1);

        c1 = _mm_loadu_ps(Ctemp+i4);
        i_col2 = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j1));
        _mm_storeu_ps(Ctemp+i4,c1);

        c1 = _mm_loadu_ps(Ctemp+i8);
        i_col3 = _mm_loadu_ps(A+((i8)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col3,j1));
        _mm_storeu_ps(Ctemp+i8,c1);

        c1 = _mm_loadu_ps(Ctemp+i12);
        i_col4 = _mm_loadu_ps(A+((i12)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col4,j1));
        _mm_storeu_ps(Ctemp+i12,c1);

        

        //unroll j
	c1 = _mm_loadu_ps(Ctemp2+i);
        //i_col = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j2));
        _mm_storeu_ps(Ctemp2+i,c1);

        c1 = _mm_loadu_ps(Ctemp2+i4);
        //i_col = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j2));
        _mm_storeu_ps(Ctemp2+i4,c1);

        c1 = _mm_loadu_ps(Ctemp2+i8);
        //i_col = _mm_loadu_ps(A+((i8)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col3,j2));
        _mm_storeu_ps(Ctemp2+i8,c1);

        c1 = _mm_loadu_ps(Ctemp2+i12);
        //i_col = _mm_loadu_ps(A+((i12)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col4,j2));
        _mm_storeu_ps(Ctemp2+i12,c1);

         

	c1 = _mm_loadu_ps(Ctemp3+i);
        //i_col = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j3));
        _mm_storeu_ps(Ctemp3+i,c1);

        c1 = _mm_loadu_ps(Ctemp3+i4);
        //i_col = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j3));
        _mm_storeu_ps(Ctemp3+i4,c1);

        c1 = _mm_loadu_ps(Ctemp3+i8);
        //i_col = _mm_loadu_ps(A+((i8)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col3,j3));
        _mm_storeu_ps(Ctemp3+i8,c1);

        c1 = _mm_loadu_ps(Ctemp3+i12);
        //i_col = _mm_loadu_ps(A+((i12)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col4,j3));
        _mm_storeu_ps(Ctemp3+i12,c1);

        

	c1 = _mm_loadu_ps(Ctemp4+i);
        //i_col = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j4));
        _mm_storeu_ps(Ctemp4+i,c1);

        c1 = _mm_loadu_ps(Ctemp4+i4);
        //i_col = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j4));
        _mm_storeu_ps(Ctemp4+i4,c1);

        c1 = _mm_loadu_ps(Ctemp4+i8);
        //i_col = _mm_loadu_ps(A+((i8)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col3,j4));
        _mm_storeu_ps(Ctemp4+i8,c1);

        c1 = _mm_loadu_ps(Ctemp4+i12);
        //i_col = _mm_loadu_ps(A+((i12)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col4,j4));
        _mm_storeu_ps(Ctemp4+i12,c1);

        

      }
    }
  }
} else {
  blocks8(m,jboun,kboun,iboun,A+ii+kk*m,A+jj+kk*m,C+ii+jj*m);
}
}
}
}

}


   void single_edgecase(int m, int jboun, int kboun, int iboun, float *A, float *A2, float *C)
  {
     for(int k = 0; k<kboun; k++)
       for (int j = 0; j<jboun; j++)
         {
         for( int i = 0; i < iboun; i++ )
            {
            C[i+j*m] += A[i+k*m] * A2[j+k*m];
            }
         }
  }

 void blocks8(int m, int jboun, int kboun, int iboun, float *A, float *A2, float *C) {
  __m128 j1,j2,j3,j4,c1,i_col1, i_col2, i_col3, i_col4, i_col5;

  int i4;
    for(int kk = 0; kk < kboun; kk+=8)
       for(int jj = 0; jj < jboun; jj+=8)
           for(int ii = 0; ii < iboun; ii+=8) {
		int iboun8 = min(8, (iboun-ii));
 		int jboun8 = min(8, (jboun-jj));
 		int kboun8 = min(8, (kboun-kk));
		if(kboun8==8&&jboun8==8&&iboun8==8) {
			for( int k = kk; k < kk+8; k++ )
  			{
    
    		for( int j = jj; j < jj+8; j+=4 )
   		 {
   		   int stored = j+1, stored2=j+2, stored3=j+3;
   		   j1 = _mm_load1_ps(A2+(j+k*m));
   		   j2 = _mm_load1_ps(A2+((stored)+k*m));
   		   j3 = _mm_load1_ps(A2+((stored2)+k*m));
  		    j4 = _mm_load1_ps(A2+((stored3)+k*m));

		   float *Ctemp  = C+(j*m);
                   float *Ctemp2 = C+(stored*m);
                   float *Ctemp3 = C+(stored2*m);
                   float *Ctemp4 = C+(stored3*m);

      
      			for( int i = ii; i < ii+8; i+=8 )
      			{
        		//C[i+j*m] += A[i+k*m] * A[j+k*m];
			i4=i+4;
        		c1 = _mm_loadu_ps(Ctemp+i);
        i_col1 = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j1));
        _mm_storeu_ps(Ctemp+i,c1);

        c1 = _mm_loadu_ps(Ctemp+i4);
        i_col2 = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j1));
        _mm_storeu_ps(Ctemp+i4,c1);


        

        //unroll j
	c1 = _mm_loadu_ps(Ctemp2+i);
        //i_col = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j2));
        _mm_storeu_ps(Ctemp2+i,c1);

        c1 = _mm_loadu_ps(Ctemp2+i4);
        //i_col = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j2));
        _mm_storeu_ps(Ctemp2+i4,c1);


         

	c1 = _mm_loadu_ps(Ctemp3+i);
        //i_col = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j3));
        _mm_storeu_ps(Ctemp3+i,c1);

        c1 = _mm_loadu_ps(Ctemp3+i4);
        //i_col = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j3));
        _mm_storeu_ps(Ctemp3+i4,c1);


        

	c1 = _mm_loadu_ps(Ctemp4+i);
        //i_col = _mm_loadu_ps(A+(i+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col1,j4));
        _mm_storeu_ps(Ctemp4+i,c1);

        c1 = _mm_loadu_ps(Ctemp4+i4);
        //i_col = _mm_loadu_ps(A+((i4)+k*m));
        c1 = _mm_add_ps(c1, _mm_mul_ps(i_col2,j4));
        _mm_storeu_ps(Ctemp4+i4,c1);


        

      }
    }
  }
		} else {
		  single_edgecase(m,jboun8,kboun8,iboun8,A+ii+kk*m,A2+jj+kk*m,C+ii+jj*m);
 		}
	   }
 }
    


